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Abstract 


This research is involved with the implementations of advanced computational schemes 
based on large eddy simulations (LES) and direct numerical simulations (DNS) to study the 
phenomenon of mixing and its coupling with chemical reactions in compressible turbulent 
flows. In the efforts related to LES, we have initiated a research program to extend the 
present capabilities of this methods for the treatment of chemically reacting flows, whereas 
in the DNS efforts, we have focused on detail investigations of the effects of compressibility, 
heat release and non-equilibrium kinetics modeling in high speed reacting flows. Our 
efforts to date, have been primarily focused on the simulations of simple flows, namely 
homogeneous compressible flows and temporally developing high speed mixing layers. 

This report provides a summary of our accomplishments during the second six months of 
the research supported under Grant NAG 1-1122 sponsored by NASA Langley Research 
Center administrated by Dr. J. Phil Drummond of the Theoretical Flow Physics Branch. 
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1. Summary of Accomplishments to Date 


During the second half of the first year of activities under this research, we have been 
primarily involved with the continuations of our tasks in the same two directions as that 
indicated in our first report. Namely: (1) development of subgrid closures for LES of 
compressible reacting flows, and (2) Understanding the mechanisms of mixing and reaction 
in high speed combustion. In the efforts related to the first task, our efforts are concentrated 
on the simulations of homogeneous reacting turbulence, whereas in the second task direct 
numerical simulations of a high speed reacting mixing layer is the subject of main focus. 
Below, a summary of our efforts during the second 6 months of this research is provided 
and detail descriptions are furnished in Appendix I and Appendix II. 

1.1. Large Eddy Simulations of Compressible Reacting Flows: 

Our major goal in this effort is to initiate a program to extend the present capabilities 
of LES for the treatment of chemically reacting flows. In the efforts to date, we have 
been primarily concerned with a priori analysis of subgrid fluctuation in a compressible 
homogeneous flows. This analysis is mainly involved with constructing the shapes of the 
simulated PDF’s within the subgrid. This work is currently in progress, and Appendix 
I and Appendix II provides a write-up on some of our recent progress. Our conclusions, 
some of which are indicated in these appendices, are summarized below: 

The primary objective of the research conducted to date has been to develop and evaluate 
mathematical models based on the PDF methods for predicting the phenomena of mixing 
and chemical reaction in turbulent reacting flows. The research was initiated in two parallel 
but related directions: 

1. Development of improved turbulence closures based on PDF methods for the purpose 
of predicting the statistical behavior of reacting turbulent flows, and also for statistical 
treatment of the behavior at the subgrid. 

2. Implementation of Direct Numerical Simulations (DNS) to provide data for evaluation of 
the performance of the closures for both Reynolds averaging and subgrid closures. 

Based on our efforts thus far, we have been able to draw the following conclusions: 

1. In the context of single-point PDF formulation, the mapping closure of Kraichnan (1989) 
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provides the most satisfactory results. 

2. The main deficiency of the above mentioned mapping closure (or any other single-point 
PDF model) is associated with its incapability in predicting the frequency (or length) scale 
of PDF evolution. 

3. The EDQNM spectral closure provides a reasonable means of overcoming the deficiency 
mentioned above in transport of conserved scalar quantities. 

4. The extension of the mapping closure for PDF formulation at two-point is very difficult. 
Among the other currently available alternatives, the LMSE closure is the simplest to use 
but does have drawbacks. 

5. Preliminary results obtained by means of two-point PDF formulation based on the LMSE 
model and the EDQNM closure indicate that the evolution of the relevant length scales is 
somewhat insensitive to the rate of chemical reaction in a chemistry of the type A + B — ► 
Products. Therefore, the EDQNM of a conserved scalar can be used in conjunction with 
the mapping closure at single-point for both non-reacting and reacting scalars. 

6. With the utilization of the mapping closure, the rate of fuel consumption in plug flow 
reactors can be predicted by a closed form algebraic relation. The results obtained by this 
relation compares very well with DNS data for a stoichiometric mixture under equilibrium 
condition. In fact, to the best of our knowledge, this relation is more satisfactory than any 
other suggested correlations (empirical or theoretical) available in the literature within the 
past thirty years! 

7. The extension of the model to account for non-stoichiometric mixtures has been done. It is 
shown that for such mixtures, the consumption rate of fuel can not predicted by a simple 
algebraic relation; rather, it can only be expressed in terms of complex integrals. 

In addition to the primary conclusions listed above, we have also been able to show that: 

8. The /3-density does reasonably well in approximating the single point PDF distribution of 
a conserved scalar property, if the evolution of its variance is known a priori. 

9. With the implementation of the /3-density, it is also possible to develop a closed form 
algebraic solution for the rate of reactant conversion under the assumptions indicated in 
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item (6). 


10. The solution corresponding to that in item (9) for non-stoichiometric mixtures would also 
be in the form of complex integrals. 

11. Despite its simplicity, the /3-density is only applicable for description of “simple” distribu- 
tions. It is not applicable for bimodal distributions, and it can not predict the subsequent 
evolution of the PDF, if the distribution is not of /3-type initially. Furthermore, it can not 
be extended for predictions of non-equilibrium reactions and there are no firm mathemat- 
ical proofs to justify its use. 

1.2. Direct Numerical Simulations of High Speed Reacting Mixing Lavers: 

This work is involved with the assessment of the roles of compressibility and heat release, 
and the non-equilibrium effects of chemical kinetics on the compositional structure of the 
flame in high speed reacting mixing layers. This work has been mainly motivated due to 
the interest of NASA in understanding the global and detail mechanisms of mixing and 
chemical reactions in high speed flows. 

Our current efforts in this work are directed toward understanding one important phe- 
nomenon. Namely, the exact influences of the Damkohler number on the statistical com- 
positional structure in a compressible mixing layer under the influence of non-equilibrium 
exothermic chemical reactions. It must be indicated that our previous work on examining 
the role of Damkohler number in incompressible flows showed that the effects of this pa- 
rameter is rather insensitive to the procedure by which it is varied. Our current efforts, 
however, indicates that this is not the case in compressible reacting layers. We do not 
have a write-up of our results addressing this issue. Nevertheless, Appendix III provides a 
summary of some of our related work in this area. 
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2. Published Work 


The following publications and reports have resulted from our efforts during the second 
six months of this research: 

2.1. Conference Papers: 

C. J. Steinberger, “Model Free Simulations of a High Speed Reacting Mixing Layer,” 
Presented at the AIAA Northeastern Regional Student Conference, held at Worcester 
Polytechnic Institute, Worcester, MA. April 15, 1991. 

2.2. Honors: 


Mr. Craig Steinberger was the winner of the first award in Graduate Full Paper Compe- 
tition at the AIAA Northeastern Regional Student Conference, held at Worcester Poly- 
technic Institute in Worcester, MA (April 1991). He has been invited to participate at 
the National AIAA Student Conference to be held at Reno, Nevada during the AIAA 
Aerospace Sciences Meeting in January 1992. 
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Appendix I 


This Appendix is provided by Dr. Cyrus K. Madnia, a Post-Doctoral Research Associate 
working on this project. 

SUMMARY 

Direct Numerical Simulations (DNS) are performed to investigate the phenomena of mix- 
ing in a decaying two-dimensional, homogeneous turbulent flow under non-reacting and 
reacting unpremixed conditions. In the reacting case, a second order stoichiometric irre- 
versible reaction of the type A + B — ► Products is considered. An infinitely fast chemical 
reaction rate is assumed without including the effects of non-equilibrium chemistry. The 
data obtained by DNS is utilized to examine the behavior of the scalar field in a stochas- 
tic manner. This involves the extraction of Probability Density Functions (PDF’s) of the 
scalar field from the DNS results. The examination of the data in this manner provides a 
useful mechanism for assessing the range of validity of PDF methods for turbulence mod- 
eling in Reynolds averaged equations, and for subgrid closures in Large Eddy Simulations 
(LES). 

Simulations are performed of flows with various levels of compressibility to examine its 
effects on the structure of the flow field. In the low compressibility regime, the flow 
exhibits features similar to those observed in previous incompressible simulations, whereas 
in the high compressibility case, shocklets are formed. In both cases, the simulated results 
indicate that the PDF of a conserved Shvab-Zeldovich scalar quantity, characterizing the 
mixing process, evolves from an initial double-delta distribution to an asymptotic shape 
that can be approximated by a Gaussian distribution. During this evolution, the PDF 
can be approximately described by a Beta density and also exhibits features compatible 
with those predicted by the mapping closure of Kraichnan. The main difference exhibited 
between the two cases is the time scale of decay of the rms of the scalar by which the PDF 
is characterized. A discussion of the results is provided with a focus on their relevance 
towards turbulence modeling and LES closure. A brief review is also provided of recent 
related contributions in PDF modeling and LES in order to provide a framework for 
anticipated continuation of the present work. 
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1. INTRODUCTION 


In review of recent computational work on turbulent reacting flows, Drummond (1991), 
Oran and Boris (1987) and Givi (1989) indicate the need for further development in both 
the methodology and the implementation of Direct Numerical Simulations (DNS) and 
Large Eddy Simulations (LES). According to these reviews, the results of previous and 
ongoing works towards the development of advanced numerical algorithms for simulating 
reacting turbulent fields have been successful, thus warranting their continued utilization 
for future implementations. With the broad knowledge gained to date, it is now widely 
accepted that foreseeable developments in advanced computational facilities will not be 
sufficient to relax the restriction of DNS to flows having small to moderate variations 
of the characteristic length and time scales (Rogallo and Moin, 1981; Reynolds, 1990; 
Schumann and Friedrich, 1986, 1987, Givi, 1989). Nevertheless, DNS can be (and have 
been extensively) used to enhance our understanding of the physics of chemically reacting 
turbulent flows by providing specific information concerning the detailed structure of the 
flow, and by furnishing a quantitative basis for evaluating the performance of turbulence 
closures. In this sense, DNS have provided an effective tool for such studies and have 
established useful guidelines for future investigations. 

LES appear to provide a good alternative to DNS for computing flows having ranges of 
parameters similar to those encountered in practical systems (Reynolds, 1990; Hussaini 
et al ., 1990, Ferziger, 1981, 1983, Schumann and Friedrich, 1987). This approach has 
a particular advantage over analogous Reynolds-averaging procedures in that only the 
effects of small-scale turbulence have to be modelled. However, despite their success for 
the analysis of both incompressible (Ferziger, 1981) and compressible (Erlebacher et al., 
1987) flows, they have seldom been employed for calculations of reacting turbulence. This 
is primarily due to appearance of enormous unclosed terms in the transport equation 
governing the evolution of the filtered quantities in LES. In these equations, the correlations 
amongst the scalar variables must be accurately modelled in order to provide a realistic 
estimate of the reactant’s conversion rates. 

A methodology which has proven useful for predictions of turbulent reacting flows is based 
on the Probability Density Function (PDF) or the joint PDF of the scalar variables which 
characterize the field (O’Brien, 1980; Pope, 1979, 1985, 1990a). A particular advantage of 
the method is the fundamental property of the probability density, namely that it provides 
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the complete statistical information about the behavior of the variable. This, in comparison 
with moment methods, provides a unique advantage in that once the distribution of the 
PDF is known, the mean and higher order moments of the reaction rate (or any other 
functions of the scalar field) can be directly determined without a need for further closure 
assumptions. The implementation of PDF methods in conjunction with the Reynolds- 
and/or Favre-averaged form of the equations of turbulent combustion has proven very 
successful (For an excellent recent review see Pope (1990a)). During the past 15 years, 
these methods have been used for simulating a variety of reacting flow configurations and, 
in many cases, the predicted results have been indicative of the higher accuracy of these 
methods in comparison with conventional moment methods. Based on this indication, it is 
safe to state that presently, from the statistical point of view, PDF methods axe the most 
attractive choice for the description of chemically reacting turbulent flows. 

A similar scenario may be surmised to be the case when the PDF methods axe used in the 
context of subgrid closures for LES. In this case, the PDF (or the joint PDF) of the scalar 
variables within the ensemble of grids (identifying the large scales of the flow) may uniquely 
determine the subgrid average value of the reaction rate (or other functions of the scalar 
variables). This PDF, in conjunction with an appropriate model for the hydrodynamic 
fluctuations within the subgrid, may facilitate the implementation of LES for the analysis 
of reacting turbulence. 

Our main goal in this work is to initiate an assessment of the use of PDF methods for their 
possible future implementation in LES, as well as to further our investigation of the PDF 
methods in modeling turbulent combustion problems. The approach followed is similar 
to those previously used in the a priori analyses. Namely, the data obtained by DNS are 
utilized to construct the appropriate PDF’s that can be used for either turbulence modeling 
or subgrid closures. We are currently in a preliminary stage of our work, especially in tasks 
directed at subgrid modeling. 

In order to establish a framework for the discussion of the results, a review is provided 
of PDF methods and subgrid closures in Section 2. This review is rather brief and is 
presented solely for the purpose of providing a reference for the remainder of the chapter. 
However, the reader is referred to other collective works and survey articles which discuss 
these statistical methods in greater detail. Also, the discussion on PDF methods focuses on 
the treatment of scalar quantities alone, without including the treatment of scalar-velocity 
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PDF’s. The specification of the problem under current scrutiny is provided in Section 3, 
in which the relevance of some of the characteristic variables are expressed, together with 
a description of the geometrical configuration and the initialization procedure. The results 
of simulations are presented in Section 4, in which we place an emphasis on the use of 
stochastic methods for modeling of the mixing phenomenon. This chapter ends in Section 
5, where we summarize our findings, and also provide some speculations and suggestions 
for future work, some of which are currently under way. 

2. BACKGROUND REVIEW 

2.1. PDF METHODS: 

In the majority of previous work on the computational treatment of turbulent combustion, 
statistical methods have generally been employed. These methods involve the represen- 
tation of physical variables in a stochastic sense, and in conjunction with appropriate 
transport equations, predict the approximate behavior of the flow field. In the framework 
of such methods, the physical quantities are treated as random variables, and the trans- 
port equations describing the evolution of these quantities are used in such a way as to 
yield their stochastic variations. At present, two general types of strategies are identified 
which can be utilized for statistical treatment of turbulent combustion problems (Libby 
and Williams, 1980): (1) the moment method, and (2) the Probability Density Function 
(PDF) method. In the first approach, “averages” (or “means”) of the relevant physical 
variables are introduced, and the governing transport equations are subsequently averaged. 
These averages are generally defined in such a way to yield the desired mathematical prop- 
erties and are usually introduced in the form of time-, ensemble-, or Favre-averages. The 
second approach, the PDF approach, relies directly on the PDF of these variables by which 
the desired moments can be obtained. 

Moment methods usually involve the determination of the statistical means by virtue of 
solving transport equations for the various unclosed moments of the equations. These 
equations provide, in a sense, a description of the variables at one level higher than the 
equations for the mean values (first moments). This results in explicit transport equations 
for the second order moments. These equations, however, are still not in a closed form, due 
to the appearance of the unclosed third order correlations which require further modeling. 
Regardless of the degree of truncation of the higher moments, the closure problem always 
remains, and transport equations at any level require modeling at higher orders. 
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The primary advantage of moment methods is their relative ease of use in engineering 
applications. Most models suggested in the literature can be implemented directly into 
computer codes developed for laminar flows. This allows the engineer to economically 
obtain approximations to combustion problems for many different sets of flow conditions. 
The results obtained in this way have, in many cases, exhibited good agreement with 
experimental measurements. This has been a major factor in justifying the use of moment 
methods in turbulent combustion and it appears that these methods will remain popular 
for engineering applications for the near future. 

The closure based on the second approach, the PDF method, has proven very useful in the 
theoretical description of turbulent reactive flows (Hawthorne et ai, 1949). The philosophy 
of this approach is to consider the transport of the PDF’s rather than the finite moments 
of the statistical variables. The key advantage of the method is the fundamental property 
of the PDF, namely that the mean reaction rate (and any other functions of the scalar) 
can be evaluated directly from the PDF. Representing the scalar field involving a species 
by $ = ((f) 1 , <j > 2 , . . . ,<p' J ), and the reaction conversion rate of species a by w c ,($), the mean 
of this rate (denoted by < >) can be written as: 

/ + oo 

Vi(E)u a (~)dE (1) 

-CO 

In this equation, E represents the scalar space and V\(E) is the PDF of the scalar variable 
at a point, defined so that: 


V\ (~)d=. = Probability (,=, < 4> < r. + dE) (2) 

where dE = ( d £* , d£ 2 , . . . , dff° ) denotes an elemental hypervolume at 5. 

The scalar PDF, V\(E), defined in Eqs. (l)-(2), is in its simplest form in that it contains 
information only about the scalar variables, 4>. Also, it governs the probability distribution 
only at a single point. However, it contains much more information than is required to 
determine the mean value of any functions of the random variables 4*. Therefore, its 
evaluation in comparison to direct evaluation of < u a > (if it could have been done) 
is understandably more difficult. Nevertheless, since 'Pi(H) includes all the statistical 
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information about the scalars, its determination is in many ways more useful than that of 
the mean values. 

A systematic way of evaluating a PDF is to obtain and solve a transport equation governing 
its evolution. A transport equation for ^(E, E = £ 1 ,£ 2 ,. . . , can be constructed by 
relating the PDF to $(x,t) in terms of Dirac-delta functions (Pope, 1979, 1985, 1990a; 
O’Brien, 1980, 1981, 1986): 


=< e(x,t) > 


where g is the single-point fine grained density defined by: 


(3) 


q = n - *°) < 4 ) 

Q — 1 


From the appropriate conservation equations for the scalar field and Eqs. (l)-(4) one can 
obtain transport equation describing the evolution of the PDF. Limiting the formulation to 
that of a constant density Fickian diffusion flow, we have (see O’Brien, 1980 for a complete 
derivation): 


< Y • V S >= -JL( Ua v,) + r < > (5) 

Here, V is the velocity field and T is the diffusion coefficient of the species. From this 
equation, it can be seen that the effects of reactivity and diffusivity (the terms on the 
RHS) are to transport the PDF into the composition space (E), whereas the convective 
transport is a physical space phenomenon. Furthermore, it is shown that the influences of 
u; Q ($) in the composition space appear in closed form. However, models are needed for 
the closures of the molecular mixing term (second term on the RHS) and the turbulent 
convection term (second term on LHS). 

Since the term that includes the interaction between the scalar variables through the 
reaction rate u Q appears in a closed form in Eq. (5), in what follows we limit the discussion 
to the transport of a single reacting scalar (i.e. a = 1) which is being convected by the 
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velocity field. In this frame, the diffusive transport may be expressed as (O’Brien, 1980, 
1981): 


< >= Vx, [E(6l6)P.Ki)l (6) 

Here, E^Hi) represents the expected value of the scalar at x 2 from its known value at 
Xj. This equation explicitly reveals a fundamental difficulty associated with the PDF for- 
mulations, namely that the transport equation for a one-point PDF (V\ (£,x,t)) requires 
information about the second point, x 2 . In general, n-point PDF’s, V n (^i,^2 
Xj ,x 2 , . . . , x n ,<) require information about the n + 1th point, x n+1 . This is further demon- 
strated by considering the transport of V n in the compositional space (Ievlev, 1970, 1973; 
Kuo and O’Brien, 1981): 


where: 


C {q) = TLim x . 

-n + l 




Vx .. [^n+l|6,6,-*-,^n,X 1 ,X 2 ,...,X n )] 

-n + 1 


( 8 ) 


and E(£ n +\ |£i, • • • , in, Xj , x 2 , . . . , x n ) is the expected value of £ n +i at the point x n+1 

based on all the values of £ 1 ,^ 2 » • • • ? in at x l5 x 2 , . . . ,x n . The dependence of C ^ on the 
n + 1th point in the transport equation for V n indicates the need for a closure model which 
relates V n +i to V n . 

Despite the attractiveness of the PDF methods in obtaining a closed representation of the 
reaction term, there is another difficulty with the implementation of the method. This 
difficulty is caused by the increase of the dimensionality of the PDF as the number of 
scalars and/or the number of statistical points increase. This imposes a severe limit on 
the maximum number n in the equation for V n . The majority of previous work on PDF’s 
have closed the equations at the first level (n = 1), and only recently has closure at the 
n = 2 level been attempted. Both these closures are reviewed in detail by Givi (1989). 


12 



2.2. LES AND SUBGRID CLOSURES: 


Despite the present capability of modern supercomputers in allowing calculations with 
more than one million grid points, the range of length and time scales that can be resolved 
by DNS is substantially smaller than those of turbulent flows of practical interest. In DNS 
the largest computable scales are limited by the size of the computational domain and the 
“turn-over” time of the large scale structures; the smallest resolvable features are limited by 
the molecular length and the time scales of the viscous dissipation and chemical reactions. 
DNS therefore, in comparison to turbulence modeling, find its greatest application in basic 
research problems in which the scales of the excited modes remain within this band of 
computationally resolved grid sizes. This implies that for an accurate simulation, the 
magnitudes of the viscosity and the diffusivity must be large enough to damp out the 
unresolved scales, and the size of the computational time step must be kept small enough 
to capture the correct temporal evolution of the turbulent flame. In this context, since the 
instantaneous behavior of the flow variables is directly simulated at all spatial points at 
every instant, the computed data can actually be used to provide a quantitative basis for 
evaluating the validity of turbulence models and for assessing the performance of subgrid 
closures in LES. 

The limitations associated with DNS may be alleviated, to an extent, by pre-filtering the 
transport equations (Aldama, 1990), a practice implemented in LES. This is effectively 
equivalent to eliminating the scales smaller than those resolvable within a given mesh. In 
this way, the variables can be represented on the number of grid points that are available 
for the simulations. This filtering procedure facilitates the simulations of flows with larger 
parameter ranges on a coarser grid. The disadvantage is that some modeling is required for 
the closure of the scales excluded by filtering. The selection between DNS and LES (and 
turbulence modeling for that matter) is dependent on the type of flow being considered, 
on the range of physical parameters that characterize the turbulent field, and the degree 
of desired accuracy. 

A straightforward implementation of LES involves the decomposition of the transport 
variables into “large scale” and “subgrid scale” components. The former is related to the 
large scale eddies in the turbulent field, whereas the latter is the component containing the 
small scale fluctuations. The pre-filtering of the variable $(x,<) is performed by means of 
the convolution integral (Ferziger, 1977, Aldama, 1990): 
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$(*,<)= IIL F£(x — x')$(x' ,t)dx' 


( 9 ) 


where F£ is an appropriate filter function with a characteristic length A, and the integra- 
tion is over the entire flowfield. The remaining portion of $ from the filtered quantity is 
defined as the subgrid-scale field and is represented by: 


$"(x,<) = $(x,<)-$(x,t) (10) 

The magnitude of pre-filtered averaged $ and its subgrid component axe obviously 
dependent on the filter type (function F£) and its size A (for reviews see Ferziger, 1981, 
1982, 1983). The simplest choice would be a box filter with: 


F£(x - x') = 


1, if |x — x'| < A; 
0, Otherwise. 


( 11 ) 


where A = (A x , A y ,A z ) represents the dimensions of the box, and the length in each 
direction should be larger than the original grid spacing in that direction. Other types of 
filters have also been suggested in the literature, each of which with computational and 
physical advantages as well as limitations (see Aldama, 1990 and Schumann and Friedrich, 
1986, 1987). 

Implementation of LES as a practical tool involves a combination of DNS for the filtered 
portion of the transport variable $, and modeling of the subgrid-scale portion, The 
transport equations for the large scale components are obtained by filtering the instan- 
taneous transport equations by means of implementing the filter defined by Eq. (9). In 
the resulting filtered equations, analogous to those in Reynolds averaged, the equations 
representing the large scale field contain unclosed terms involving the fluctuation of the 
subgrid components. The methodology practiced to date in dealing with these fluctuations 
is very similar to that employed in Reynolds averaging procedures. One may be able to 
develop and solve transport equations for the moments of fluctuations (Ferziger, 1981). 
These equations, analogously, contain some higher order moments which are needed to be 
modelled. The number of these terms are usually more than the corresponding ones that 
appear in Reynolds averaged transport, because some items that are zero in the Reynolds 
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averaging approach are non-zero when filtering is used. Correspondingly, the tasks re- 
quired for modeling the subgrid fluctuations would be somewhat more complicated than 
those in turbulence modeling. Nevertheless, since the small scales of turbulence are be- 
lieved to exhibit a more “universal” character is the main reason to anticipate that the 
approach based on subgrid modeling would be more promising than the procedures based 
on ensemble averaging closures. 

Within the past two decades, there have been many attempts to fine-tune the models by 
optimizing the closure for subgrid fluctuations and the type of filtering to identify the large 
scale components of the variables. In most of these efforts, the closure of hydrodynamic 
fluctuations has been the subject of major concern. There has also been some results 
for the passive scalar simulation and modeling of velocity-scalar correlations. However, 
no attempts have been made for the treatment of reacting scalars and the treatment of 
scalar-scalar interactions. This is probably due to the consensus among the researchers in 
this field that until an accurate subgrid model is constructed to represent the evolution of 
non-reacting scalar variables, the extension to reactive flow simulations will be a difficult 
task. 

With the new developments and progress in PDF methods and the advantages offered 
by such schemes over moment methods, one may anticipate that these methods may also 
prove useful in LES. In this case, the influences of the chemical reactions and the subgrid 
scalar-scalar correlations can be included by means of considering the PDF’s of the fluc- 
tuating scalars within the computational grid. The advantage of using these methods for 
the subgrid closure is apparent since once these PDF’s are known, any statistical quan- 
tity related to scalar fluctuations can be subsequently determined. This determination, 
although an ambitious task, is not impossible. The obvious and simplest choice is to follow 
an approach based on guessed PDF methods. Similar to Reynolds averaging, the first two 
subgrid moments of the scalar variables can be solved by LES and then the shape of the 
PDF can be parameterized based on these two moments. This parameterized approxi- 
mate distribution can be appraised by a comparison with the “exact” shape constructed 
via DNS results. Similar to common practice in a priori assessments (Erlebacher et a/., 
1987), the PDF distribution can be specified by performing two sets of calculations; one 
with a coarse mesh, the other with a fine mesh. The results obtained from the fine mesh 
simulations can be used to construct the PDF distribution which in turn may be used in 
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extensive subsequent simulations on the coarse mesh. In this setting, calculations with 
large transport parameters which otherwise could not be simulated on the coarse grids, 
are possible. 

A more direct, but substantially more complicated procedure involves the solution of trans- 
port equation for the PDF’s of the subgrid scalars rather than assuming their form. This 
approach, like its counterpart in turbulence modeling, has the advantage that the effect 
of scalar-scalar correlations appear in a closed form. However, models are needed for the 
molecular diffusion within the subgrid, and one has to resort to mixing models (a subject 
of current intense research) for providing the closure. The approach based on guessed 
PDF methods is feasible and is within the reach of present day computers. The approach 
based on solving the transport equation for the subgrid scale PDF might require extensive 
computational resources (Pope, 1990a). These models must initially be developed in a 
simple flow. A homogeneous flow is an excellent configuration for this purpose, and will 
be discussed in the next section. After the establishment of a successful model, it may be 
utilized for simulating more complicated flows. The extent of this utilization is dependent 
on the available computational resources and on the performance of the model in rigorous 
trials. 


3. DESCRIPTION OF THE PROBLEM 

The subject of our DNS is a two-dimensional homogeneous box flow under the influence 
of a binary chemical reaction of the type A + B — * Products. To impose homogeneity, 
periodic boundary conditions are employed in all directions of the flow, identifying the box 
as a homogeneous member of a turbulent universe. This periodicity allows the mapping 
of all the aerothermochemical variables from the physical domain into a Fourier domain, 
thus allowing the implementation of spectral methods for numerical simulations. The flow 
field is assumed to be homogeneous and isotropic, and is initialized in a similar manner 
to that of Passot and Pouquet (1987). This involves the superposition of a “random” 
velocity to a zero mean velocity and also includes random initial density, temperature, and 
Mach number fluctuations. These fluctuations have certain energy spectra, and the ratio 
between the compressible and the incompressible kinetic energy can be varied to assess 
the effects of compressibility. The specification of the fluctuating field with a random field 
is to mimic a “probabilistic” turbulent field in the context of “deterministic” DNS. This 
approach is similar to that followed in previous direct simulations of turbulent reacting 
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flows (Riley et al . , 1986; Givi and McMurtry, 1988; McMurtry and Givi, 1989). In a 
laboratory flow, these fluctuations appear as the result of interactions with the surrounding 
universe. Such interactions do not exist in the isolated homogeneous flow considered in 
our simulations. Therefore, in order to introduce the “noise,” which plays a central role 
in laboratory turbulence, these perturbations are randomly superimposed to initialize the 
fluctuating field for DNS. The generated turbulence field is of decaying nature, i.e. there 
is no artificial forcing mechanism to feed energy to the small wave numbers. While it 
is more desirable to have a stationary turbulence field to focus only on the behavior of 
the scalar field (Eswaran and Pope, 1988; McMurtry and Givi, 1989), it is not clear how 
to implement such a forcing mechanism in the compressible case without modifying the 
behavior at low wave numbers. 

The scalar fields are defined to be square waves with the two species out of phase and 
at stoichiometric conditions. The species field is assumed to be dynamically passive, in 
that turbulence influences the consequent transport of the scalar field with the neglect of 
reverse influences. In the non-reacting case, the trace of only one of these reactants is 
considered, whereas in the reacting case, the transport of an appropriate Shvab-Zeldovich 
variable is assumed to portray the reactant’s conversion rate. This is possible by assuming 
an infinitely fast chemical reaction rate, and by assuming that the reactants have identical 
thermochemical properties. In this framework, the effects of non-equilibrium chemistry 
are neglected; the inclusion of such effects are postponed for future investigations. 

The computational package employed in simulations is based on the modification of a 
code developed by Erlebacher et al. (1987). This code is based on a spectral collocation 
algorithm with Fourier trial functions. All the variables are spectrally approximated on 
N 2 collocation points, where N represents the number of collocation points in each of 
the directions. The spatially homogeneous flow evolves in time, and at each time step 
N 2 defines the sample data size for the purpose of statistical analysis. A third order 
accurate Runge-Kutta finite difference scheme is employed for temporal discretization. 
For a trustworthy simulation, the magnitudes of the Reynolds and Peclet numbers must 
be kept at moderate levels and the size of the time step should be kept small. The code is 
capable of simulating both two- and three- dimensional flows and we have performed both 
such simulations. In the presentation of our results in the next section, however, we limit 
our discussion to that of a two-dimensional flow. 
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4. PRESENTATION OF RESULTS 


Computations are performed on a domain with a normalized dimension of 2ir in each of 
the directions of the flow. With the available computational resources, a resolution of 
256 x 256 collocation grid is attainable. With this resolution, the magnitude of the Taylor 
microscale Reynolds numbers that could be simulated is in the range, .Re* « 20 ~ 30. 
Simulations are performed with a wide spectrum of compressibility levels; here we only 
report the results obtained by use of two extreme cases: one with a low compressibility 
level, or pseudo incompressible, and the other with a relatively high level. In the former, 
the initial value of the normalized density rms is very small, i.e. p r ma = y / < P 12 >/po & 0, 
whereas in the latter, this value is of order unity. The compressibility level was monitored 
by adjusting the initial values of the following parameters (Passot and Pouquet, 1987): 
(1) the rms of the Mach number, Mt = v< M' 2 >, and (2) the ratio of the compressible 
energy to that of the total kinetic energy, Below, we limit the discussion to the results 
obtained for two cases. In what follows, we refer to pseudo incompressible simulations in 
which the following values are used initially: M t = 0.2, x — 0.01, and to compressible 
simulations in which at the initial time Mt = 0.6, x = 0.2. 

The compressibility effects can be manifested by both flow visualization and by considering 
the global behavior in an integral sense. In the former, the contour plots of the relevant 
variables constructed from the DNS results show the qualitative behavior, whereas in the 
latter, the ensemble averages of the simulated results portray the quantitative response. 
To demonstrate this, in Figs. 1 and 2 we present the contour plots of the density for 
the pseudo incompressible case and for the compressible case, respectively. A prominent 
difference between the two cases is the formation of steep gradients in the high compressible 
case which are not observed in the pseudo incompressible simulations. The regions of high 
gradients are referred to as shocklets, and consistent with the previous simulations of Passot 
and Pouquet (1987), appear when the initial levels of density and Mach number fluctuations 
are high. The effects of increased compressibility can be quantitatively demonstrated by 
examining the temporal variations of the maxima and minima of the divergence of the 
velocity, the fluctuating Mach number and the ratio of the compressible to the total kinetic 
energy. These are presented in Figs. 3-6. In Figs. 3 and 4 the temporal variations of the 
minimum and maximum values of the divergence of the instantaneous velocity (V • Y) are 
presented for the pseudo incompressible and for the compressible cases, respectively. Note 
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that in the low compressible flow, the minimum and maximum values are approximately 
mirror image of each other (with respect to V • Y = 0). In the high compressible case, 
however, the maximum and minimum divergence curves are asymmetric. The minimum 
value of divergence occurs at the normalized time of f ~ 0.7. This time corresponds to 
the onset of formation of the shocklets, and indicates severe compression within the flow. 
It must be mentioned that the capture of these shocklets by means of global numerical 
methods (without any artificial dissipation) is rather difficult. With 256 x 256 collocation 
points, this was the strongest shock that we were able to capture, and a lower resolution 
would result in a significant oscillation in the magnitude of the velocity divergence after the 
appearance of the shock. One must be careful not to confuse these numerical oscillations 
with compression and expansion waves. 

The effects of compressibility are further quantified in Figs. 5-6. For the pseudo incompress- 
ible flow, Fig. 5, the magnitude of M t decreases monotonically, whereas in the compressible 
case, Fig. 6, the decrease in M t is interrupted by plateaux. The locations of these plateaux 
coincide with those corresponding to a rise in the compressible kinetic energy and the local 
minima of the velocity divergence. Therefore, these locations correspond to the times at 
which shocklets can be formed in the flow. In the pseudo incompressible flow, the mag- 
nitude of x remains fairly constant, indicating that the initial low level of compressibility 
remains low at all times. 

With the development of the flow, the species field would consequently evolve. To visualize 
the flow, the contour plots of species .4 are presented in Fig. 7. Parts (a) and (b) of this 
figure correspond to the zero and the infinitely fast reaction rate cases, respectively. This 
figure exhibits the effects of random motion on the distortion of the scalar field and the 
mixing of the two initially segregated reactants (the contours form parallel lines at t = 0). 
The effects of chemical reaction are to increase the steepness of the scalar gradients and 
to reduce the instantaneous values of the reactants, as indicated by a comparison between 
parts (a) and (b) of the figure. The quantitative behavior of the scalar development is 
depicted in a statistical sense by means of examining the evolution of the PDF’s of the 
conserved Shvab-Zelclovich variable J . This variable is normalized and defined within the 
region [0,1]. Correspondingly, its PDF, T\(J) is always bounded in this region. The 
temporal variation of V\{J) is presented in Fig. 8. It is shown in this figure that at 
the initial time, the PDF is approximately composed of two delta functions at J = 0,1, 
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indicating the two initially segregated reactants A and B. At later times, it evolves through 
an inverse-like diffusion in the composition space. The heights of the delta functions 
decrease and the PDF is redistributed at other J values in the range [0, 1]. At even later 
times, the PDF becomes concentrated around the mean value. Proceeding further in time 
results in a sharper peak at this mean concentration, and the PDF can be approximated 
by a Gaussian distribution. This Gaussian type behavior has been observed in previous 
simulations of Givi and McMurtry (1988), and Eswaran and Pope (1988), and also has 
been corroborated in a number of experimental investigations. 

An interesting character of these PDF’s is that throughout their evolution, the simulated 
results compare remarkably well with that of a Beta density. This is also shown in Fig. 8 by 
a comparison between the Beta density and the DNS generated PDF’s. The Beta density 
is parameterized with the same first two moments as those of the DNS. Therefore, in all the 
figures, the results are presented with respect to a time scale (f*) proportional to the decay 
of the variance of the scalar J . Higher order moments of the DNS data for the variable J 
also show good agreements in comparison with those predicted by the Beta density. This 
is demonstrated in Figs. 9 and 10, where the normalized kurtosis (/x 4 ) and superskewness 
(^g ) are presented of the random variable J and the reactant A, respectively. At time 
zero, these moments are close to unity and monotonically increase as mixing proceeds. 
For the Shvab-Zeldovich variable J, the magnitudes of the kurtosis and superskewness 
resulting from the Beta density asymptotically approach the limiting values of 3 and 15, 
respectively. These correspond to the normalized fourth and sixth moments of a Gaussian 
distribution as the variance of tends to zero (as t* ► oo). The DNS generated results 
are very close to those of the Beta distribution throughout the simulations. However, the 
limiting value for variance of J approaching zero can not be obtained in the simulations 
due to obvious numerical difficulties. In the reacting case, the moments of the scalar A are 
consistently higher than those of the conserved scalar, but portray similar trends in both 
Beta density and DNS generated results. 

The trends shown above are also observed in the compressible simulations. The profiles of 
the PDF’s and those of the higher order moments and their comparisons with the corre- 
sponding parameterized Beta density are presented in Figs. 11-13. Again, the comparison 
is remarkably good. The main difference between these results and those of pseudo incom- 
pressible simulations is the time scale of the decay of the variance by which the PDF is 
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evolved, f This time scale can not be incorporated into the PDF description in the format 
employed here. The single-point nature of these PDF’s do not allow for any information 
about the length scales (or any other scales requiring two-point statistics) of turbulence. 
For a systematic inclusion of the length scale into the PDF and quantitative description 
of the differences between the results in the two cases, one is required to consider the evo- 
lution of the PDF’s at two points (at least). Employing the results of DNS to evaluate (or 
to generate) models for more than one-point-level closures (O’Brien, 1985; Eswaran and 
O’Brien, 1989) has proven to be a challenging task, and is the subject of current research 
(Jiang and Givi, 1991). 

The results of these simulations indicate that the approximation of a Gaussian PDF for 
the final stages of mixing of a conserved scalar is well justified. This corroborates the 
results of laboratory experiments (e.g. Miyawaki et al., 1974; Tavoularis and Corrsin, 
1981) and those of other numerical simulations (e.g. Eswaran and Pope, 1988; Givi and 
McMurtry, 1988). The numerical experiments performed here, however, axe similar to most 
of the laboratory investigations in that they include the results of only one experimental 
realization. Future numerical experiments using different initializations and/or with better 
numerical resolution and including three dimensional effects would be logical extensions of 
this work. 

In addition to the Beta distribution approximations, the DNS results also compare very 
well with the mapping closure recently introduced by Kraichnan (see Chen et al., 1989 
and Kraichnan, 1990). This closure has the property of allowing the relaxation of the 
PDF of a conserved scalar property to an asymptotic Gaussian distribution. Pope (1990b) 
has utilized this closure for the prediction of PDF evolution in an isotropic homogeneous 
turbulent flow with an initial double delta distribution (similar to the condition imposed 
here). He showed that in the context of a one point PDF, the results obtained by this 
closure compare favorably with the DNS results. However, again there is no length scale 
information built into the model, and only the evolution of the PDF without any informa- 
tion regarding the time scale of evolution can be predicted. Madnia and Givi (1991) have 
extended this model for predicting the decay rate of a reacting scalar in homogeneous tur- 
bulence. The generated results obtained by applying this model for the prediction of the 

f This is not illustrated in the figures clearly, since the time is scaled by the rate of variance 
decay. 
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decay rate of a reacting scalar for the cases considered in DNS are shown in Figs. 14 and 15. 
In these figures, the DNS generated results and those predicted by a Beta density are also 
shown. The two PDF’s (by the mapping closure, and the Beta density) are constructed in 
such a way to yield the same first two normalized moments of the Shvab-Zeldovich variable 
as those of DNS. In this context, the predicted results of the decay by both models show 
good agreements. The agreement for the mapping closure model is rather surprising, con- 
sidering that this closure was developed for an isotropic, fully three-dimensional turbulent 
flow. This is probably due to matching of the normalized first two moments. A better 
test of the model would be its utilization at two-point level, and its comparison with DNS 
generated results. Also, the agreement between the mapping closure results and those of 
DNS worsens as the compressibility level is increased. This again is understandable in that 
Kraichnan’s model was developed for a purely incompressible flow without including any 
compressibility effects. 

With the construction of the DNS database, the results are used to construct the PDF’s 
of the scalar quantities within the subgrid. Since the major riddle in PDF modeling 
is associated with the closure of the diffusion term and not the chemical reaction term, 
we limit the analysis to that of the generated PDF’s of the conserved Shvab-Zeldovich 
variable. The construction of the PDF’s is implemented by placing a coarse (n 2 ) mesh 
over the physical space occupied by the fine mesh (IV 2 ). This mesh can be considered as 
the one to be potentially used in LES. In this way, it is assumed that the LES predicts the 
averaged results at the centers of the n 2 mesh, and the fluctuations within the coarse grids 
are to be modelled by an appropriate PDF closure. The actual LES in this case would 
correspond effectively to simulations with a homogeneous box filter in which the values of 
the filtered means are constant. This is somewhat analogous to the procedure followed by 
Schumann (1989). The main difference is the mechanism by which the fluctuations are to 
be considered. Schumann (1989) simply neglected such effects, but with consideration of 
the PDF, it may be possible to account for such effects, albeit statistically. Within each 
cell of the coarse mesh, we have the values of the scalar quantities at N# = ( N/n ) 2 equally 
spaced points. This means that there is an ensemble of N<j, sample points at each location 
to construct the PDF’s. This construction is implemented by statistical sampling of N# 
number of data points. The domain of the scalar property 4> € [<j> min ,4> max] is divided 
into a number of bins of size A <f>, and the number of scalar variables within each of these 
bins are counted. The PDF of the variable <j>, denoted by V(() is calculated in the interval 
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£ < <t> < ( + d£ = d<j>, from the definition of the probability: 


V(£)d£ = Probability (£ < 4> < £ + d£) 


( 2 ') 


which translates into: 

V(£)d£ = Probability Density 

Number of elements in the interval £ < <f> < f + d( (12) 

= Ni 

The value of N+ is a measure of the width of the filter, in that it determines the scale 
at which the fluctuations of the scalar field are not accounted for deterministically, but 
are included in a probabilistic manner by the PDF’s. Consequently, the shapes of these 
distributions are dependent on the magnitude of the sample size N If this value is too 
small, then there will be a large scatter in the data. If it is too large, then the PDF’s 
correspond to those appropriate for Reynolds averaging and not LES. In this primeval 
implementation, we have not yet performed an actual LES (using a PDF subgrid closure). 
Rather, we have focused only on constructing the shapes of the PDF’s from the DNS 
results. These PDF’s can again be quantified by the magnitudes of their moments to 
provide a reasonable a priori assessment of the closure for future actual LES. 

With the 256 2 available total data, there is a limit on the magnitude of the ensemble data 
for an a priori assessment. Our experience indicated that an ensemble of about 32 2 was 
the lowest acceptable limit for a meaningful statistical analysis. Therefore, the domain was 
broken into 8x8 squares, within each the PDF’s were constructed. There was also some 
analysis with 4x4 squares for better statistics. However, since this is very close to the 
Reynolds averaging limit, the generated results will not be considered in the forthcoming 
discussions. 

The PDF’s constructed from this ensemble showed that as the mixing proceeds, all the 
PDF’s tend to have an asymptotic Gaussian distribution. This is to be expected, since 
a subset constructed from a large random set with a probability distribution tends to 
portray the same statistical behavior. However, the results showed that for all the cases 
considered, the PDF’s tend to have an evolution similar to that of the Beta density. This is 
demonstrated in a sample plot of the distributions in Fig. 16 for both pseudo incompressible 
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and compressible simulations. In this figure, the PDF’s within one of the coarse grids are 
compared with that of a Beta density with the same two first moments as those within the 
subgrid. The comparison is fairly good, considering the size of sampled data. The overall 
deviation from the Beta density and the Gaussian distribution in all the subgrids can be 
quantified by measuring the norms of the deviation of higher moments. The L 2 norms of the 
fourth moments for the pseudo incompressible and the compressible simulations are shown 
in Figs. 17 and 18, respectively. Note that at the initial time, the deviation from a Beta 
distribution is very small and increases with progress in time, but not substantially. At 
final times, the DNS data are closer to that of a Gaussian distribution, but the magnitudes 
of the second moments are not small enough for the Beta density to approach a Gaussian 
distribution. Similar overall trends are also observed for the compressible simulations. 

With the results of the analysis, albeit primitive, it is speculated that the approach based 
on PDF parameterization from its first two moments may prove serviceable for LES, at 
least for a conserved scalar. However, the approach requires the accurate input of the 
first two moments of the subgrid which must be provided by modelled LES transport 
equations. In the absence of a better alternative, the Beta density may prove useful, and 
the generated results can be used to estimate the maximum reaction conversion rate in the 
infinitely fast chemistry limit. However, the generalization to include finite rate kinetics 
may not be very straightforward, since the first two moments (including the covariance of 
the scalar fluctuations) must be known (or modelled) a priori. By the same token, the 
extension to more complex kinetics with the inclusion of non-equilibrium effects, which axe 
very important for practical applications, remains to be a challenging task. 

An improvement of the procedure obviously involves the solution of a transport equation 
for the PDF (rather than its parameterization with its moments). However, the com- 
putational burden may be prohibitive due to the increased dimensionality of the subgrid 
PDF transport equations. An estimate of the computational requirements indicates that 
the cost associated with the implementation of LES procedure involving the stochastic 
solution of the PDF transport equation may be of the same order as that of DNS on the 
fine grids, unless the ratio of the fine to coarse grid is large. In this regard, the tradeoff 
between LES and DNS depends on other factors, such as the physical complexity and 
computational resources. 
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5. CONCLUDING REMARKS 


A spectral collocation algorithm has been employed to study the mechanism of mix- 
ing and reaction in a non-premixed, decaying two-dimensional homogeneous turbulent 
flow. The evolution of the species field in a one-step stoichiometric reaction of the type 
A + B —> Products was the subject of the investigation. Calculations were performed for 
zero-rate and infinitely fast rate chemistry in both pseudo incompressible and relatively 
high compressible flows. The effects of compressibility are delineated by the formation of 
thin shocklets and can be quantified by a global examination of the data. Mixing and the 
reaction conversion rate (in the reacting case) are characterized by examining the evolu- 
tion of a Shvab-Zeldovich conserved scalar quantity extracted from the DNS. The results 
show that the PDF of this scalar variable evolves from an initial double-delta function 
distribution to an asymptotic shape which can be approximated by a Gaussian distribu- 
tion. During this evolution, a Beta density qualitatively describes the DNS generated 
PDF’s in both pseudo incompressible and compressible simulations. This is quantitatively 
demonstrated by a good agreement between the values of the higher order moments of the 
PDF’s calculated from the DNS data and those predicted by the Beta approximation. The 
generated results also compare very well with those predicted with the mapping closure 
of Kraichnan. The decay of the mean of the reacting scalar determined by DNS compares 
very well with that predicted by this model and that obtained via the use of the Beta den- 
sity. This trend is observed in both pseudo incompressible and compressible simulations. 
The main difference is the time scale of the decay of the variance of the PDF’s. 

The Beta density seems to also predict the behavior of the fluctuating field within the 
subgrid for a conserved scalar quantity. Therefore, it may provide a reasonable means of 
parameterizing the subgrid fluctuations in a stochastic sense. The procedure, however, 
requires the input of the first two moments of the variable within the subgrid, and these 
must be provided by modelled LES transport equations. Also, the generalization to include 
finite rate kinetics may not be very straightforward, since these moments must be modelled 
a priori. An improvement of the procedure involving the solution of a transport equation 
for the PDF is possible, but may not be practical due to the increased dimensionality of 
the modelled subgrid PDF transport equations. From the computational point of view, 
it is recommended to assess the applicability of such a procedure in (yet) simpler flows, 
before its implementation in simulating more complex flows. 
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In comparing the DNS results with those of predictions, the first two moments of the Shvab- 
Zeldovich variable must be matched. This is because single-point PDF’s do not include 
information about the frequency scales of turbulence. Therefore, the evolution of the length 
scale, or any other parameter requiring two-point information, must be introduced in an 
ad hoc manner (here, such information is provided via DNS). Future validation studies of 
the PDF’s by DNS need to consider the evolution of multi- (at least two-) point statistics. 
O’Brien (1985) has already employed a two-point formulation in the statistical description 
of homogeneous turbulent flows. The preliminary results obtained by such formulations 
are in good agreement with experimental data (Eswaran and O’Brien, 1989). However, a 
comparison between the predicted results and the DNS data, similar to the ones reported 
here for single point PDF’s is recommended. 

Future work should also deal with finite rate chemistry with the inclusion of non equilibrium 
effects, and also on examining the extension of PDF methods in dealing with the subgrid 
closures in LES of such flows. 
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Figure 1. Plot of density contours for the 
pseudo incompressible case, t = 6.142. 



Figure 2. Plot of density contours for the 
compressible case, t = 0.686. 





Figure 4. Temporal variation of the minimum and maximum values of 
the velocity divergence for the compressible case. 




t 


Figure 5. Temporal evolution of M t (solid line) and 
X (dashed line) for the pseudo incompressible case. 
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Figure 6. Temporal evolution of M t (solid line) and 
X (dashed line) for the compressible case. 







Figure 7. Plot of species A concentration contours at 
t* = 1.0715. (a) non-reacting case; (b) reacting case. 



Figure 8. Temporal evolution of V\(J) for the pseudo incompressible 
case. DNS data (dotted line), Beta density (solid line). 

(a) f = 0.035, (b) = 0.549, (c) t* = 0.800, 



































Figure 12. Temporal variation of the kurtosis, ^ 4 , for 
the compressible case, (a) variable J , 

(b) Species .4 in reacting case. 









Figure 13. Temporal variation of the superskewness, /X6, for 
the compressible case, (a) variable J , 

(b) Species A in reacting case. 
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Figure 14. Mean concentration of the reacting species A 
vs. time for the pseudo incompressible case. 
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Figure 15. Mean concentration of the reacting species A 
vs. time for the compressible case. 



Figure 1G. Sample PDF’s of the variable J within the subgrid, 
(a) pseudo incompressible case, t * = 0.0355, (b) pseudo incompressible case, t 
(c) compressible case, t * = 0.101, (d) compressible case, t* = 1.507. 
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Figure 17. L 2 norm of deviation of the fourth moment 
vs. time for the pseudo incompressible case. 
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Figure 18. L 2 norm of deviation of the fourth moment 
vs. time for the compressible case. 




Appendix II 


The results presented in this Appendix are obtained by Mr. Alex Tsai, who is a M.S. 
candidate working on this project. 

In this Appendix, we demonstrate the capabilities of the EDQNM spectral closure in 
its implementation in a two-point PDF formulation. In doing so, we present the results 
obtained by means of the solution of a two-point PDF transport equation in a homogeneous 
turbulent flow under the influence of a binary reaction of the type A + B — * Products. The 
turbulent field is statistically stationary, with a specified energy spectrum. The EDQNM 
model is used for the closure of the convective flux, and the LMSE model is used for 
the modeling of the molecular mixing. The purpose of this appendix is to show that the 
EDQNM is indeed capable of predicting the evolution of the length scales, and that the 
chemistry does not have a significant impact on the evolution of these scales. 

In Fig. II. 1 , we present the temporal evolution of the maximum mean fuel concentration for 
two different initial values of the integral length scale. This figure shows the exact effects of 
the length scale. Note that as this scale is smaller, i.e. the reactants are in a closer vicinity 
of each other, the rate of chemical reaction become larger, and the fuel is consumed faster. 
This figure displays the “beauty” of two-point formulation! In Fig. II.2, the evolution of 
the Taylor microscale is presented for the fuel concentration under two limiting conditions 
of no chemistry and infinitely fast reaction rate. Note that until the final stages of reaction 
(at long times) the length scales are the same, and the results only deviate slightly at final 
times. The same is also evident in Fig. II. 3, in which the normalized dissipation rates of 
the scalar quantities are presented. Note that the difference between the two cases is not 
substantial. 

With the results presented above, we conclude that the evolution of the length scales 
predicted by means of the model described above, is not sensitive to the rate of chemistry. 
This is validated by considering the scales under the two limiting conditions. However, the 
model employed for the closure of molecular mixing is based on the LMSE closure. While 
the results of our preliminary DNS validate the conclusions drawn from these two-point 
simulations, we would still need to perform more simulations, in order to make this point 
more clear. 
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Figure II. 1. Temporal variation of < A >t / < A >o 
(Combination of LMSE and EDQNM Closures) 



Figure II.2. Temporal variation of the Taylor Microscale 
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Figure II. 3. Temporal variation of the Scalar Dissipation Rate 

(EDQNM Closures) 





This Appendix is prepared by Mr. Craig Steinberger, a Ph.D. candidate working on this 
projects. 
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Abstract 

The effects of compressibility, chemical reaction exothermicity and non-equilibrium 
chemical modeling in a combusting plane mixing layer were investigated by means of 
two-dimensional model free numerical simulations. It was shown that increased com- 
pressibility generally had a stabilizing effect, resulting in reduced mixing and chemical 
reaction conversion rate. The appearance of “eddy shocklets” in the flow was observed 
at high convective Mach numbers. Reaction exothermicity was found to enhance mix- 
ing at the initial stages of the layer’s growth, but had a stabilizing effect at later times. 
Calculations were performed for a constant rate chemical rate kinetics model and an 
Arrhenius type kinetics prototype. The Arrhenius model was found to cause a greater 
temperature increase due to reaction than the constant kinetics model. This had the 
same stabilizing effect as increasing the exothermicity of the reaction. Localized flame 
quenching was also observed when the Zeldovich number was relatively large. 


1 Introduction 

Research on the subject of compressible reacting mixing layers has been of high priority in 
recent years. Much of this effort has been devoted to the development of high speed air 
breathing flight vehicles. This type of vehicle would, according to current proposals, use 
supersonic combustion ramjet (scramjet) engines for propulsion. In such an engine, fuel 
is injected into a high speed airflow. The mechanisms of mixing and combustion of this 
non-premixed, high speed, compressible, reacting flow is of great complexity. 

The purpose of this study was to investigate the coupling between the hydrodynamic 
and chemical phenomena that are present in compressible, reacting mixing layers. Although 
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experimental studies have provided useful results, this avenue of research is often difficult 
and expensive. Numerical simulation has become widely accepted as an alternative and 
supplement to experimental investigation. Model free simulations have proven especially 
useful because there is no need to introduce the potential inaccuracies of a turbulence model. 
An extensive review of current numerical simulation techniques and results for reacting flows 
is provided by Givi [1], 

1.1 Previous Research 

Brown and Roshko [2] found that the turbulent mixing layer is dominated by large scale 
coherent structures, or vortices. These structures convect at a nearly constant speed and 
tend to coalesce with neighboring vortices. Papamoschou and Roshko [3] continued these 
experiments to determine the effect of compressibility on the spreading rate of a supersonic 
mixing layer. They found that it is useful to study the flow in reference frame that travels 
with the flow at the same speed as an average large scale structure. A parameter which 
quantifies the compressibility in the flow was proposed as the convective Mach number, 
Af c ; defined as the Mach number of the flow with respect to the above mentioned frame of 
reference. A direct correlation was found between M c and the stability of the flow. Lele [4] 
replicated the results of [3] by means of direct numerical simulation of a two-dimensional 
mixing layer. In addition, he noted the appearance of eddy shocklets within the flow for 
M c > 0.7. 

McMurtry, et al. [5] studied the effects of an exothermic chemical reaction on the shear 
layer via direct numerical simulation. In this study, an approximate set of equations that 
are valid for low Mach number flows was used. It was found that the heat liberated from 
a chemical reaction causes the layer to grow at a slower rate that of on-heat releasing flow. 
These results agree with those obtained experimentally. In addition, direct numerical sim- 
ulations concerned with compressibility and heat release effects on a compressible mixing 
layer were performed by Givi et al. [6]. 

The phenomenon of flame extinction in non-premixed flames has been the topic of theo- 
retical and experimental study, although direct numerical simulations of such flows have been 
somewhat limited. Recent reviews of some of the prevalent theories regarding the structure 
of turbulent non-premixed flames, as well as some of the experimental and numerical work 
in that field are provided by Bilger [7] and Peters [8]. 



1.2 Scope of Present Research 


The scope of the present work is to examine the effects of compressibility and chemical 
reaction exothermicity on a reacting plane mixing layer. An examination is also made of 
the non-equilibrium effects of the chemical kinetics on the structure of a flame. These are 
accomplished by direct numerical simulation of an unsteady two-dimensional shear layer. 
The governing equations are integrated via high order finite difference methods, and no 
turbulence modeling is employed. Since the focus of this study is the general physics of this 
type of flow, physical modeling is kept as simple as possible so that the above mentioned 
hydrodynamic and chemical effects may be isolated. The use of simplifying assumptions has 
the added advantage of saving considerable amounts of computational resources, without 
reducing the accuracy or validity of the results. 

2 Method of Investigation 

A two-dimensional version of the SPARK computer code [9, 10] was used to investigate 
a chemically reacting, compressible, temporally developing mixing layer. The full Navier- 
Stokes equations as well as the applicable chemical species conservation equations (see [9] 
or [11] for more detail) were solved via numerical integration. The methods used were the 
Gottlieb- Turkel finite difference scheme [12] and a dissipative compact parameter method 

[13]. Both of these algorithms are of fourth order spatial accuracy and second order temporal 
accuracy. 

One of the primary assumptions made in these simulations is that the mixing laver is 
temporally developing. That is, the reference frame of the simulations is defined to be moving 
with a velocity equal to that of the mean flow. The advantages of this approximation arl 
twofold. First, with the temporal assumption the inflow and outflow boundary conditions can 
be assumed to be periodic. The eliminates the problem of specifying downstream boundary 
conditions. Second, the temporal assumption means that only a relatively small region of the 
flow is being simulated. This region is then followed in a Lagrangian sense as time p-ogresses. 
This results in considerable savings in computational resources, which means that the flow 
may be simulated in greater detail. 

The chemical reaction in the flow is assumed to be of a simple, irreversible, second order 
type of the form 

A + B — ► Product -f Heat nj 

The reaction is characterized by the kinetics mechanism, which is given by the single step 



model of 

w-KjCaCb ( 2 ) 

where Ca and Cb represent the concentrations of the reacting species and are assumed equal 
at the free streams, i.e. Ca^ — Cb x = C x . Kj is the reaction rate constant, and can be 
normalized to form the definition of the Damkohler number, Da: 


Da = — ^ Cco 

f/oo / 6 W |o 

In the present study two types of chemistry models 
constant Kj) and an Arrhenius type model in which 
is written as 


( 3 ) 

were used; constant rate kinetics (i.e. 
Kj varies with the temperature. This 


Kj = Ajt °° 


( 4 ) 


where Aj is the pre-exponential factor and Ze is the Zeldovich number, defined as 


Ze = 


E 

RToo 


( 5 ) 


Here E is the activation energy and R is the universal gas constant. is the free stream 
temperature. When the Arrhenius kinetic model is used, the pre-exponential factor Aj 
replaces Kj in the definition of Da (Eq. 3). 

Combustion exothermicity is measured by the energy liberated by the chemical reaction, 
AH°. The magnitude of this energy is parameterized by a non-dimensional heat release 
parameter Ce, defined by: 


Ce = 


-AH 0 

CvToo 


( 6 ) 


Thus, Ce = 0 corresponds to a non-heat releasing chemical reaction. 

The flow was initialized with a hyperbolic tangent streamwise velocity profile, with re- 
actant A on the top half of the layer with a free stream velocity of +f/co and reactant B 
on the bottom half with a free stream velocity of — This flow configuration is shown 
schematically in Fig. 1. There was no initial fluid motion in the transverse direction, and 
the pressure was assumed to be initially constant throughout the flowfield. Depending on 
the problem simulated, the temperature was either assumed to be initially constant or had 
an initial Gaussian distribution. In most cases, the mixing layer was perturbed by numerical 
truncation errors. However, when necessary, harmonic forcing was explicitly added to trigger 
the flow instabilities. For a majority of the simulations, a 128 x 128 point grid was used, 
with grid compression in the areas of the flow with the sharpest gradients. In those cases 



where the gradients were too strong for the details of the flow to be resolved, a 256 x 256 
point grid was used. 

All parameters are normalized when appropriate by initial or free stream conditions. 
Time is normalized by 

1 = Woo (7) 


3 Sample Results 

3.1 Effects of Compressibility 

The effect of compressibility on the mixing layer was investigated by varying M c , the convec- 
tive Mach number over a range of values from M c = 0.2 to M c = 1.2. The vorticity thickness 
vs. normalized time for different values of M c is given in Fig. 2. This figure clearly shows that 
the growth rate of the vorticity thickness is decreased with increased compressibility. Com- 
pressibility also affects the time needed for the layer to roll up into a vortex, which is shown 
in the figure as a jump in vorticity thickness. Valuable information may be obtained from 
examination of the variation in the transverse direction of the normalized average streamwise 
component of the velocity and the normalized mean square streamwise velocity fluctuations. 
The average streamwise velocity profile is shown in Fig. 3. The most significant feature of 
these curves is the steepness of the mean velocity profiles at high convective Mach numbers. 
This shows a sharper velocity gradient across the layer, implying a lesser rate of mixing. 
The degradation of mixing has the effect of retarding the progress of the chemical reaction. 
The suppression of turbulence fluctuations with compressibility is shown by the profiles of 
the mean square of the fluctuating velocity in Fig. 4. A marked decrease in the amplitude 
of the fluctuations can be observed. 

Another feature of the increased compressibility is the formation of eddy shocklets as a 
result of increased compression after the formation of large scale structures. These shocklets 
are formed for both subsonic and supersonic free stream flows in order to adapt to high 
pressure at the braids. An example of this can be seen upon examination of Mach number 
contours, shown in Fig. 5. This phenomenon has not been observed in any of the experiments 
on compressible shear layers, to date. This might be caused by the fact that the simulations 
are two-dimensional, whereas a laboratory flow is dominated by small scale three dimensional 
transport. It is anticipated that future three dimensional simulations will clarify the matter. 



3.2 Effects of Reaction Exothermicity 

Calculations were also performed to assess how reaction exothermicity affects the mixing 
layei. Results are presented of simulations for constant rate kinetics with heat release values 
of C e — 0, 1.5, and 6. One set of calculations was performed using Arrhenius type kinetics 
with Ce = 1.5 and Ze — 10. 

The results of these simulations presented in the form of plots of the voriticitv thickness 
versus normalized time (Fig. 6). This figure shows that the rate of growth is highest for the 
no-heat release case (C e = 0), and, as the heat release parameter is increased, the growth 
of the layer slows. The relatively smooth regions of voriticity thickness growth mav be 
attributed to diffusion thickening, and a jump in vorticity thickness represents vortex roll- 
up. For the C e = 0 case, the layer responds to background perturbation; fairly quickly, and 
vortical structures are formed at t m ss 3. An increase in the magnitude of the heat release 
results in a delay of vortex roll-up, and the jump in vorticity thickness does not occur until 
t ~ 7. Further increase of the heat release parameter results in additional delays, as can 
be seen for the case of C e = 6. This behavior is also observed for the Arrhenius model with 
Ce — 1-5. In these two cases, the effects of exothermicity are most pronounced; vorticity 
roll-up does not occur at all, and the only growth in the thickness of the mixing layer is due 
to molecular diffusion. 

The influence of the heat release on the structure of the flame is demonstrated by examin- 
ing the product thickness of the layer, defined by the normalized total product concentration 
of the layer as a function of time (Fig. 7). At initial times, the effect of heat release i; an 
enhanced product formation, while the reverse applies at the intermediate and final stages. 
Initially, the effect of heat release is to expand the fluid at the cores of the laver. A large 
mixing zone is formed, which results in an enhanced reaction and an increased product 
formation. This explains the increased initial product thickness. However, as the heat re- 
lease increases and the layer thickens, the growth of the instability modes become subdued, 
postponing the formation of the large scale vortices. After initial times, the non-heat release 
(C e = 0) and C e = 1.5 cases predict a sharp increase in the product thickr.ess. Thi; is caused 
by the dynamics of the large scale structure formation. The mechanism i: roll-up causes the 
leactants to be entrained into the layer from their prospective free streams. Thi; produces 
an increase in the extent of mixing, reaction, and product formation. The lack of roll-up in 
the C e =6 and Arrhenius cases means the only mixing process is molecular diffusion. Since 
in these simulations the diffusion mechanism is not as efficient as the roil-up of vortices in 
the mixing of reactants, the product thickness is correspondingly lower. 



3.3 Flame Extinction 


Another area of study was the nonequilibrium effects leading to local flame extinction. Sim- 
ulations were run for large values of the Zeldovich number for an Arrhenius kinetics model 
and for a constant rate kinetics chemistry model as a control. To obtain the vortex roll-up 
necessary to this analysis, explicit harmonic forcing was added to the initial conditions of 
the flow. A forcing amplitude of 0.5% of the mean velocity was used. The most unstable 
frequencies were found by a linear stability analysis of an equivalent incompressible flow. It 
was found that if the flow was excited only by the most unstable frequency, two vortices 
form but they do not pair. If the first subharmonic is added, the vortices undergo pairing. 

In the case of large Zeldovich numbers, it was observed from the simulations that, after 
roll-up and pairing, the flame undergoes local extinction at the braids. This behavior is 
depicted in contour plots of the reaction rate at times t* = 1, 1.5, 2.25, and 2.5. These 
contours are shown in Fig. S for time t’ = 1, when the vortices have just formed. The reaction 
rate is highest along the mixing surface of the layer, that is, in the center of the intervortex 
braids. The contours at time t~ = 1.5 (Fig. 9) portray the behavior at the initial stages ofS 
pairing. The reaction rate has begun to decrease in the braids of the emerged vortex. At 
time t‘ = 2.25, shown in Fig. 10, the layer has completed the pairing process and strong 
gradients are observed in the vortex braids. Since the braids are “stretched” as the vortex 
rotates, they are the area of the highest strain. The flame at the onset of extinction is shown 
in Fig. 11. Note that the reaction rate has gone to zero at the braids. Product concentration 
contours (not shown) indicate that very little product exists in these extinguished regions, 
demonstrating that the flame did not quench due to depletion of reactants. At this point in 
time, the flame is not continuous, forming what will be called for lack of a better term, a 
“flame eddy.” 

This extinction phenomenon has been explained by Peters [8j as follows: At the regions 
of high strain, the reactants are supplied at a faster rate than they can be consumed by 
the flame. Thus the local temperature in that area drops and the flame becomes very 
rich with reactants. As a result, the flame is quenched in that area. If a fast chemistry 
model or an equilibrium chemistry model is used, the extinction mechanism can not be 
captured. Investigation of such phenomena requires finite-rate chemistry simulation such as 
that presented here. 



4 Conclusions 


Direct numerical simulations of a compressible, temporally developing, reacting plane mixing 
layer were performed. Several simplifying assumptions were made so that the effects of the 
variation of isolated parameters could be studied in detail. In particular, the chemical reac- 
tion was assumed to be of the general form, A + B — » Products + Heat, and thermodynamic 
properties were assumed constant and identical for all species. 

Studies of various flow phenomena were performed by varying one representative nondi- 
mensionalized parameter while keeping all other parameters constant. The convective Mach 
number, M c , was used to describe compressibility, the heat release parameter, C e , repre- 
sented the exothermicity of the reaction, and the Damkohler number and the Zeldovich 
number described the rate of the reaction. 

The simulations concerning the effect of compressibility on the mixing layer showed a 
direct correlation between increased compressibility and increased stability (along with re- 
duced turbulence). When the convective Mach number was increased, the rate of growth of 
the mixing layer, which was measured by the growth of the vorticity thickness, was markedly 
reduced. Degradation of the development of the streamwise fluctuating velocity and mean ve- 
locity profiles was also noted. For high compressibility cases, “eddy shocklets” were observed 
within the flow. This has been previously reported in two-dimensional simulations, but has 
not been observed in experiments or in three-dimensional studies to date. An expansion of 
this work into three dimensions is suggested for future work in this area. 

Increased exothermicity was observed to slow the growth of large scale structures. At the 
initial stages of development, high heat release increased the amount of product formed via 
volumetric expansion of the core of the layer. However, the heat release caused the layer to 
be Aless responsive to perturbations, reducing the growth of the layer at later stages. The 
overall effect of increased heat release, therefore, was to stabilize the flow and to decrease 
the extent and efficiency of the reaction. 

The selection of chemical kinetics model was shown to have significant effect on the 
development of the flow. The introduction of an Arrhenius chemistry model had a stabilizing 
effect, thus degrading the progress of the reaction. Non-equilibrium chemistry effects due 
to the chemistry model was also studied. It was found that at high Zeldovich numbers, the 
flame would be quenched at regions with large local values of the strain rate. 
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Figure 1: Schematic diagram of a temporally evolving mixing layer. 



Figure 2: Normalized vorticity thickness versus normalized 
convective Mach number. 
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Figure 6: Normalized vorticity thickness versus normalized time for different values of the 
heat release parameter. 











